The user who contributed the test problem is mentioned by name. I assume it's ok if you have their permission?

Would you have time to chat about the test? I don't suspect a problem, but it would help me understand some of the details of the setup, and it would also be good to have a look at the results

A nice cross-check might be to instead use ESL to add/remove heat to/from the model at the same total rate as in the current setup and make sure you get the same results for temperature. (Addition/removal of an amount of heat from the aquifer should be the same as production/decay of the same amount of heat distributed between the water and solid because the water and solid equilibrate anyway)

Another could be, starting with the existing setup, add/remove all the heat to/from just the water, and in another run all to/from just the solid, to make sure you get the same temperature results (assuming the same total rate of heat addition/removal in each case)

Is it worth clarifying that the decay in each phase is a rate per amount of that phase (not per volume of aquifer), or is that obvious enough? A similar question applies to straight-up "decay" (sorption aside) in GWT, which is only in the water and is (presumably?) a rate per volume of water (not per volume of aquifer)

"The dimensions of decay for zero-order decay in the aqueous phase is energy per length cubed per time." -- Suggest modifying to "The dimensions of zero-order decay in the aqueous phase are energy per length cubed per time."
Similar suggestion for the analogous description of decay_solid.

Related question: based on the code and the test problem, the units of decay in the solid are energy per mass (of solid) per time, not "energy per length cubed per time"

Is it worth renaming subroutine est_fc_dcy to est_fc_dcy_water (similar for est_cq_dcy)?

In est_cq_dcy, remove part of comment that says solid still needs to be added

Remove temporary comment notes about first-order decay from est_fc_dcy, est_cq_dcy, and est_cq_dcy_solid. If desired, could replace them with formal comments saying that first-order decay of energy is not supported. Also looks like there's a reference to first-order decay in read_options

Remove temporary comment note about negative temps from est_fc_dcy. If desired, could replace it with a formal comment about the fact that negative temps are currently not checked for or prevented

Suggest modifying
if (.not. lname(2) .and. .not. lname(3)) then
write (errmsg, '(a)') 'Zero order decay in either the aqueous &
&or solid phase is active but zero-order rate coefficients, &
&either in the aqueous or solid phase, is not specified. &
&Either DECAY_WATER or DECAY_SOLID must be specified in the &
&griddata block.'
to this
if (.not. lname(2) .and. .not. lname(3)) then
write (errmsg, '(a)') 'Zero order decay in either the aqueous &
&or solid phase is active but the corresponding zero-order &
&rate coefficient is not specified. Either DECAY_WATER or &
&DECAY_SOLID must be specified in the griddata block.'

(shortened a bit, and corrected "coefficients ... is" to "coefficient ... is")

get_zero_order_decay appears to limit decay to as to avoid negative temps, which is should not be doing. there is also a kluge note about wanting to use different expressions for the rates. we should chat about this

I believe logical operators associate from left to right, so the following

  elseif (this%idcy == 2 .and. this%idcysrc == 1 .or. &
          this%idcysrc == 3) then ! zero order decay aqueous phase

which would evaluate the "and" first, then the "or", should be

  elseif (this%idcy == 2 .and. (this%idcysrc == 1 .or. &
          this%idcysrc == 3)) then ! zero order decay aqueous phase

Likewise,

  elseif (this%idcy == 2 .and. this%idcysrc == 2 .or. &
          this%idcysrc == 3) then ! zero order decay in the solid phase

should be

  elseif (this%idcy == 2 .and. (this%idcysrc == 2 .or. &
          this%idcysrc == 3)) then ! zero order decay in the solid phase